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THE  GRADIENT  METHOD  FOR  INTERFACE  TRACKING 


Introduction 

Interface  tracking  methods  are  used  in  numerical  calculations  to  represent  an  internal  physical 
boundary  as  a  discontinuity  that  has  important  properties  distinct  from  the  medium  on  either  side. 
To  approximate  the  interface  as  a  discontinlty,  the  physical  boundary  must  be  thin  compared  to  the 
domain  of  interest  and  the  importance  of  the  internal  structure  of  the  boundary  must  be  limited  to 
the  global  effects  which  the  interface  method  attempts  to  resolve  and  simulate.  There  are  several 
groups  of  interface  methods:  surface  tracking  methods  [1,2],  which  keep  track  of  distinct  points  on 
the  Interface,  volume  tracking  methods  [3,4,5],  which  reconstruct  the  interface  on  the  basis  of  a 
marker  quantity  assigned  to  each  computational  cell,  and  moving  grid  methods  [6,7,8],  which  alter 
the  computational  grid  to  keep  the  interface  on  cell  boundaries.  In  general,  all  of  these  methods 
break  the  interface  modeling  process  into  three  steps:  1)  defining  the  location  of  the  interface, 
2)  determining  the  movement  of  the  interface,  and  3)  incorporating  the  relevant  effects  of  the 
physical  boundary.  Useful  summaries  of  interface  tracking  are  found  in  the  reviews  by  Hyman  [9] 
and  Laskey,  et  al.  [10]  and  the  introductory  material  by  Hlrt  and  Nichols  [5]. 

This  paper  introduces  the  gradient  method,  a  new  interface  method  for  the  subset  of  interface 
problems  where  material  on  one  side  of  the  boundary  is  converted  into  the  material  on  the  other. 
The  gradient  method  differs  fundamentally  from  the  methods  previously  mentioned  in  that  it  is 
not  necessary  to  define  the  location  or  movement  of  the  interface.  Rather,  the  effects  of  interface 
propagation  are  determined  from  two  characteristics  distinct  to  the  physical  system.  First,  large 
gradients  are  present  as  a  result  of  material  conversion  at  the  physical  interface,  and  the  magnitudes 
of  these  gradients  are  used  to  identify  the  computational  cells  containing  the  reaction  front.  Second, 
the  interface  propagates  at  a  defined  speed  in  a  direction  normal  to  itself.  This  speed  is  a  global 
quantity  that  resiilts  from  the  combined  effects  of  heat  transfer,  diffusion,  chemistry,  and  other 
processes  relevant  to  the  particular  problem.  The  relatively  localized  large  gradients  are  combined 
with  the  propagation  speed  to  determine  the  amount  of  material  conversion  and  thus  simulate  the 
movement  of  the  interface.  The  effects  of  the  interface  movement  are  incorporated  without  the 
need  to  resolve  the  interface  beyond  the  accuracy  inherent  in  the  numerical  convection  algorithms. 
Below  we  show  how  the  method  is  applied  to  a  problem  with  propagating  chemical  reaction  fronts. 

Mannicripf  approved  December  3.  1987. 


The  Gradient  Method 


Consider  a  system  containing  a  mixture  of  gases,  A  and  B,  which  can  react  to  form  a  product 
species,  C, 

.4  +  S  — C.  (1) 

Reactions  can  occur  in  one  or  several  regions  in  the  system.  The  boundaries  of  these  reacting 
regions  are  reaction  fronts  with  large  gradients  in  species  concentrations.  In  the  description  that 
follows,  we  choose  the  number  density  of  product  C,  denoted  c,  as  an  indicator  of  the  reaction 
fronts.  Later  we  consider  other  possible  choices  for  the  indicator  variable. 

The  existence  of  a  reaction  front  does  not  directly  enter  into  the  governing  equations  for 
conservation  of  mass  and  momentum.  For  exothermic  reactions,  a  source  term  appears  in  the 
energy  conservation  equation,  but  this  source  term  is  a  function  of  reaction  rate  and  ultimately 
the  concentration  of  the  constituent  species.  Since  the  product  number  density  has  been  chosen 
as  the  indicator  variable,  the  conservation  equation  of  most  interest  is  that  for  c.  An  Eulerian 
representation  for  the  evolution  of  the  product  number  density  is 

dc 


oc  „ 

-  +  V.cF  =  U,e, 


(2) 


where  ^  is  the  velocity  vector  of  the  flow  field  and  is  the  production  term  due  to  the  propagat¬ 
ing  reaction  front.  The  left  hand  side  of  Eq.  (2)  can  be  solved  by  any  method  for  solving  continuity 
equations,  and  the  production  term  is  added  to  the  convection  solution  by  timestep  splitting  meth¬ 
ods.  The  gradient  method  is  used  to  evaluate  Wg.  For  simplicity,  other  effects,  such  as  molecular 
diffusion,  have  been  omitted  from  Eq.  (2),  but  these  may  also  be  included  using  timestep  splitting. 

Consider  a  long  enclosure  of  length  L  and  constant  cross  section  5  through  which  a  reaction 
front  is  propagating.  The  reactant  and  product  number  densities  of  the  enclosed  volume  vary  with 
time.  As  the  reaction  front  propagates  along  the  enclosure,  the  change  in  product  number  density, 
Ac,  is  given  by 


Ac  = 


(C2  -  Cl)  IS 
LS 


(3) 


where  ci  and  C}  are  the  local  product  number  densities  upstream  and  downstreaun,  respectively,  of 
the  reaction  front  and  /  is  the  distance  the  reaction  front  has  moved  along  the  enclosure  during  a 
given  period  of  time.  Eq.  (3)  is  valid  even  if  only  a  portion  of  the  enclosure  is  considered,  as  long 
as  that  portion  contains  the  reaction  front.  However,  L  would  then  be  defined  as  the  length  of 
the  portion  of  the  enclosure  under  consideration.  Taking  the  limit  as  the  length  L  gets  vanishingly 
small  gives 

Ac=|Vc|/,  (4) 


where  |Vc(  is  the  magnitude  of  the  gradient  of  the  product  number  density  across  the  reaction 
front.  The  distance  /  is  measured  normal  to  the  reaction  front,  which  is  the  same  direction  as  that 
of  the  gradient.  If  the  speed  of  the  front  is  the  local  burning  velocity,  then  /  may  be  approximated 
by 

/  =  ,  (5) 


where  Vj  is  the  local  burning  velocity  and  is  the  time  interval.  The  amount  of  product  formed 
locally  in  time  At  as  the  reaction  front  moves  is 

ttfcdt  =  Ac  =  V),At  |Vc| .  (6) 

The  integral  of  the  gradient  through  the  reaction  front  is  fixed  by  the  upstream  amd  downstream 
values  of  c  and  not  by  the  purely  numerical  details  of  the  gradient  resolution.  Thus,  the  integrated 
aonount  of  product  formed  is  the  same  as  that  produced  by  the  propagation  of  the  physical  interface. 
For  many  reactions  charaicterized  by  Eq.  (1)  the  burning  velocities  are  available  in  the  literature  for 
various  concentrations,  temperatures,  and  pressures.  In  a  numerical  simulation,  the  time  interval  is 
the  numerical  timestep,  which  is  determined  by  stability  criteria  of  the  numericad  method  or  other 
properties  of  the  coupled  system  of  reactive  flow  equations. 


Calculatioa'of  |  Vcj 

Let  us  denote  the  cell  centers  of  a  two-dimensional  computational  grid  as  (i,j),  where 
i  s  1, ...,n«  and  j  =  l,...,ny,  and  n,  and  riy  are  the  number  of  cells  in  the  x  and  y  direc¬ 
tions,  respectively.  Species  number  densities  are  defined  at  these  cell  centers.  A  central-difference 
approximation  to  a  numerical  derivative  on  this  grid  is  second-order  accurate  if  the  function  in 
question  varies  smoothly  over  the  three  cells  involved.  A  one-sided  difference  is  only  first-order 
accurate,  but  the  number  of  cells  over  which  c  must  vary  smoothly  is  reduced  from  three  cells  to 
two.  In  the  problem  of  a  thin  reaction  front,  we  have  found  that  a  one-sided  difference  gives  a 
better  approximation  to  the  derivative  because  the  largest  change  in  c  usually  occurs  over  only  two 
cells.  In  the  gradient  method,  both  the  left  and  right  hand  one-sided  differences  are  evaluated,  and 
the  value  chosen  for  the  x-direction  derivative  Qx  is  the  one  with  the  largest  absolute  value.  The 
same  procedure  is  used  to  determine  the  y-direction  derivative,  g^.  This  process  is  summarized: 


gxiiiJ) 


~  Cj.j 


9T2{iJ) 


~~  ^t— IJ 


(7a. 6) 
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{7c,  d) 


9yl{i,j)  = 


Vi+i  -  9} 


9y2{iJ)  = 


9]  -  9j-i 
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?xl 

li^xil  >  |yx2l 

if 

9x2 

15x2 1  >  l5xt| 

9yl 

l5yll  >  l5y2l 

if 

9y2 

l5y2l  >  l5yll 

5y(».i)  = 


Note,  Zi  and  are  the  z  and  y  locations  of  the  center  of  cell  (i,j)  ,  and  yxi,  yx2,  yyi,  and  yy2 
are  the  one-sided  derivatives  from  which  and  are  chosen.  The  gradient  magnitude  {Vc|  is  the 
magnitude  of  the  vector  sum  of  gg  and  y,, 

|Vc|=[y2  +  y2]'/^ 


•I 


ModiScatioas  to  |Vc| 

The  value  of  |Vc|  is  nonzero  whenever  the  values  of  c  in  adjacent  cells  are  not  identical.  But 
it  is  possible  for  this  quantity  to  be  nonz^o  yet  unrelated  to  the  presence  of  a  reaction  front, 
for  example,  in  a  compression  or  rarefaction  region.  Because  this  model  assumes  that  the  largest 
gradients  are  connected  with  reactions,  the  problem  of  unrelated  gradients  is  avoided  by  defining 
a  minimiiTn  gradient  for  reaction  and  ignoring  any  gradient  below  this  critical  value.  This  limit 
prevents  the  propagation  of  the  reaction  front  due  to  numerical  diffusion  and  restricts  the  thickness 
of  the  reaction  front  to  less  than  two  ceils.  Without  this  limit,  the  reaction  front  is  distributed 
across  too  wide  a  region  and  the  effective  burning  velocity  is  too  low. 

The  limiting  value  of  |Vc| ,  denoted  |Vcj^,  j  ,  is  expressed  in  the  form 

IVcUu  •  (''» 

where  Ci  is  a  constant  discussed  below  and  is  a  characteristic  length.  The  quantity  Cmax{iyj) 

is  the  maximum  possible  product  number  density  in  cell  (i,j),  and  is  calculated  as 


c(t,j)  =  Ci,j  + 


\r’  •  /».j 


where  fij  and  n.-j  are  the  fuel  and  oxidant  number  densities,  respectively,  N/  and  No  are  the 
number  of  fuel  amd  oxidant  molecules  consumed  per  reaction,  and  No  is  the  number  of  product 


molecules  produced  per  reaction.  The  reactants  A  and  B  in  Eq.  (1)  are  now  associated  with  the 
fuel  and  oxidizer. 

The  characteristic  length  l\  is  a  measure  of  the  distance  over  which  the  product  number  density 
goes  from  0  to  Cmax  when  a  reaction  front  is  present  in  cell  (i,j)  and  its  neighbors.  The  value  of  lx 
is  defined  as  the  weighted  average  of  the  distance  betxveen  the  center  of  cell  (i,j)  and  the  centers 
of  the  adjacent  cells  that  contribute  most  to  the  value  of  |Vc| , 

h{i,3)  =  {  max[(|5*i|  “  .0]  *  (li+i  -  Xj) 

+  max((|^,2|  -  l^xiD.O)  *(x,-  -  Xi_i) 

+  -  l^vzl)  .0)  *  (jj+i  -  y,) 

+  max[(|yvj|-|5yil).0)*(yy-y,_i)  }  /  {  max[(|y*il  -  |yx2l)  ,0] 

+  max[(|y*2l  -  |y*i|),0] 

+  max[(lyi,i|  -  |y„2l),0] 

+  max[(|yy2l  -  lyyiD.O]  },  (12) 

Note  that  the  denominator  in  Eq.  ( 12)  can  alternately  be  expressed  as  { |  |yxi  I  ~  \9x2 1 1  + 1  l^yi  I  ~  Ipy2 1 1 } ■ 
The  weighting  factors  in  Eq.  (12),  the  differences  between  the  one-sided  derivatives,  aire  constructed 
so  that  only  the  spacings  between  cell  centers  of  cells  containing  the  reaction  front  are  used  to 
calculate  /].  These  cells  are  identified  as  having  the  larger  of  the  magnitudes  of  9x2  a^d  the 

larger  of  the  magnitudes  of  9y\  and  gy%.  Using  the  difference  between  magnitudes  for  the  weighting 
factor  serves  two  purposes.  First,  the  difference  filters  out  any  small  uniform  gradient  present 
throughout  the  flow.  Second,  the  difference  puts  most  of  the  weight  on  the  direction  in  which  the 
gradient  is  changing  fastest.  Empirically  this  results  in  a  more  accurate  reproduction  of  the  input 
burning  velocity. 

The  second  case  in  which  |Vc|  must  be  modified  is  for  cells  in  which  reactions  must  occur 
regardless  of  the  calculated  magnitude  of  the  gradient.  For  example,  consider  Figure  1.  If  cell 
(i  —  l,j  —  1)  is  highly  reacted,  a  term  which  is  quantified  below,  then  cell  {i,j)  should  react.  Cells 
(i  -  l,j)  and  (i,J  —  1)  react  because  the  large  amount  of  product  in  cell  (i  -  l,y  -  1)  leads 
to  values  of  |Vc|  greater  than  in  these  three  cells.  But  the  amounts  of  product  in  cells 

(*  ”  l»j)  and  (i,j  —  1)  are  not  sufficient  to  produce  a  value  of  |Vc|  greater  than  i  in  cell 

(t,j).  Physically,  though,  the  reaction  should  proceed  in  the  corner  of  cell  (i,j). 

Another  case  in  which  a  cell  should  react  due  to  proximity  to  a  highly  reacted  cell  is  when 
two  reaction  fronts  are  merging.  Assume  in  Figure  2  that  cell  (i,j)  has  been  reacting  due  to  the 


5 


gradient  from  either  of  its  neighbors.  As  the  reaction  proceeds,  the  value  of  c,,j  increases  until 
the  difference  between  it  and  its  neighbors  is  not  sufficient  to  produce  a  value  of  jVcj  greater  than 
should  Continue  to  react,  though,  so  that  the  reaction  fronts  can  merge. 

It  is  for  situations  such  as  those  above  that  the  term  highly  reacted  and  a  second  critictil 
gradient  magnitude  are  defined.  This  second  critical  value,  ,  >  represents  a  second  minimum 

magnitude  of  the  gradient  for  cells  that  are  near  highly  reacted  cells.  The  measure  of  the  degree  of 
reaction  in  a  ceil  is  the  fraction  of  reaction  completed. 


Fi.i  = 


(13) 


The  maximum  value  of  F  for  the  eight  cells  neighboring  cell  (»,  j)  is  Fmaxihj)-  If  FmaxiiJ)  exceeds 
a  certain  value.  Ferity  the  cell  (i,j)  is  said  to  be  neighboring  a  highly  reacted  cell  amd  the  magnitude 
of  the  gradient  associated  with  cell  (i,j)  should  not  be  less  than  value  of  Feru  should 

not  be  less  than  Ci  or  else  ail  cells  would  react  before  |Vc|  exceeds  i .  The  value  of  , 

is  defined  as 


cjij) 


(14) 


where' C]  is  a  constant  and  /]  is  a  characteristic  length.  The  length  /}  is  similar  to  li  except  that 
the  weighting  factor  here  is  the  amount  that  F  for  a' neighbor  cell  exceeds  the  value  Ferit-  This 
change  in  weighting  factor  is  because  |Ve|^,  j  is  looking  for  high  gradients  to  indicate  reaction, 
but  |Ve|^.f  j  is  responding  to  highly  reacted  cells. 

The  characteristic  length  fj  depends  on  whether  the  highly  reacted  cell  is  adjacent  to  a  side  of 
cell  (t,  j)  or  on  a  diagonal.  For  the  adjacent  cells. 


=  {  max^fi+xj  -  iVtt) ,0]  ♦  (*i+i  -  Xi) 

+  max[(fi_ij  -  F’cr<«).0]*(xi  -x<_i) 

+  max[(Fi.j+x  -  i’crit)  ,0]  •  (yy+i  -  yj) 

+  max((/i,j_i  -  Ferit), 0]  * (y>  -  y,_x)  }  /  {  max[(Fi+ij  -  Ferit)  ,0] 

+  max((Fi_i  , j  ~  F ^it ) ,  0] 

+  max[(F,i+i  -  Ferit)  ,0] 

+  max[(Fi.>_x  -  Ferit),0]  }  (15) 


Here  a  characteristic  length  l-i.adj  corresponds  to  the  maximum  F  value  of  the  adjacent  cells.  .A, 
similar  calculation  relates  a  characteristic  length  h.diag  a-nd  a  corresponding  maximum  F  to  the 
diagonal  cells. 


~  {  niaX  -fcrtt)  »  0]  * 

+  max[(f’i+i,j+i  -  /’crtt),0)  ♦  dsij 
+  max[(/’i+i,j-i  -  Fcrit)  ,0)  *  dsij.i 
+  max  [(/i-ij+i  -  Ferit) .  0)  *  dsi^ij  }  / 

{  max[(/'i_i,j_i  -  f’cr<t),0] 

+  max[(f’i+i,j+i  -  ,0] 

+  max[(/’j+ij_i  -  Fcrtt)  ,0] 

+  max[(fi_i,j+i  -  /’crtt)  ,0]  } 

and 


(17) 


^ma*,diag  “  max(/' i_i  j— i  ,F Fi^ij—iy  F i+l,j+-l)  • 


Here 


=  [(*<+1  -  *.)*  +  (yy-t-i  -  yj)*] 


1/3 


s  distance  between  centers  of  diagonal  neighbors. 


(18) 


(19) 


The  characteristic  length  I2  is  either  l2,a4i  h,diag,  depending  on  whether  Fmax.adj  or  F-nax.diag 
exceeds  the  value  Fcrit-  If  both  Fmax,adi  and  Fmax,dU9  exceed  f'crtt.  then  hihj)  is  set  to  l2,adj 
because  an  adjacent  cell  will  have  greater  influence  on  its  neighbor  tham  a  diagonal  cell. 


{h,adii*tj)  F inax,adj  >  ^ rrit 

if 

h,diag{*ij) 


■  max,diag  ^  "  crit 


>  Ferit  and  F, 


max,adi  ^  FcrU 


(20) 


and 


Fmax  (itj)  ~  OlAX(Fmax,adj  »  F naXidiagi  • 


(21) 


The  discontinuity  of  Eq.  (21)  when  Fmax,adj  =  Fmax.diag  is  not  important  because  this  condition 
does  not  arise  in  practice.  Also,  tests  show  that,  except  during  early  stages  of  a  calculation, 
I7  —  h.adj  whenever  |Vc(^^j  2  is  used.  Finally,  if  cell  {i,j)  is  highly  reacted  but  all  of  its  neighbors 


are  not,  then  the  resulting  values  from  Eqs.  (15)  and  (17)  are  zero.  To  insure  that  Eq.  ( 14)  yields 
a  usable  value,  is  restricted  to  a  value  at  least  as  large  as  the  smcdlest  spacing  to  a  lujighbor 

cell. 


h(iyj)  =  niajc{(2  from  eq.  (20), 

min[(i<+i  -i,),(x.-  -Xi_i),(yj+i  -  yj)  ,(yj  -  yj-i)]}-  (22) 


Modification  for  noaunifocm  temperature  field 


The  propagation  of  an  unconiined  flame  front  is  approximately  a  constant  pressure  process, 
but  the  release  of  energy  introduces  large  temperature  gradients  at  the  front.  If  the  process  is 
to  remain,  isobaric,  the  temperature  change  must  be  reflected  in  a  proportional  change  in  species 
densities.  This  process  occurs  in  addition  to  changes  in  species  density  due  to  chemical  reactions. 

The  present  model  assumes  that  large  differences  in  the  product  number  density  are  the  result 
of  chemical  reactions  and  the  amount  of  reaction  is  proportional  to  the  difference  in  product  number 
densities  in  adjacent  cells.  To  avoid  having  density  differences  reflecting  temperature  differences 
rather  than  reactions,  all  calculations  are  done  using  number  densities  scaled  to  an  arbitrary  ref¬ 
erence  temperature,  TVe/.  This  scaling  is  done  assuming  that  all  species  obey  the  ideal  gas  law. 
Thus,  the  scaled  number  density,  is  given  by 
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The  final  value  of  ^  is  scaled  back  to  the  local  cell  conditions  by  multiplying  by  the  factor 
Trt/ITij.  It  may  be  possible  to  avoid  this  modification  if  the  model  is  formulated  in  terms  of 
gradients  of  mole  fraction  rather  than  gradients  of  number  density. 


Summary  of  the  calculation  of  |Vc| 

The  calculation  of  |  Vcj  follows  the  following  steps: 

1.  Scale  the  product  number  densities  to  the  reference  temperature. 

2.  Calculate  ( Vc| . 

a.  Calculate  the  sample  gradients  in  the  x-direction,  gsiii,j)  and 

b.  Set  the  x-direction  gradient  for  a  cell,  gs{i,j),  equal  to  the  sample  x-direction  gradient 
with  the  largest  absolute  value. 

c.  Calculate  the  sample  gradients  in  the  y-direction,  9y2{iij)- 


d.  Set  the  y-direction  gradient  for  a  cell,  gy{i,j),  equal  to  the  sample  y-direction  gradient 
with  the  lairgest  absolute  value. 

e.  Calculate  the  magnitude  of  the  gradient,  |  Vc| ,  equal  to  the  vector  sum  of  its  components. 

3.  Determine  the  minimum  gradient  magnitude  for  reaction,  j  . 

a.  Calculate  maximum  possible  product  number  density  in  a  cell,  Cmax{i,j)- 

b.  Calculate  the  characteristic  length  l\. 

c.  Calculate  the  value  of  i 

4.  Determine  the  maximum  value  of  the  fraction  of  reaction  completed,  Fmaxihj)  for  sa^ch  cell. 

a.  Calculate  the  fraction  of  reaction  completed  for  each  cell,  Fij. 

b.  Determine  the  maximum  value  of  F  for  cells  adjacent  to  cell  (i,  j). 

c.  Determine  Frnax  ,diag(iij)i  the  maximum  value  of  F  for  cells  on  the  diagonal  of  cell 

d.  Assign  the  largest  of  Fnax,adj(iJ),  Fmax,diag(iJ),  and  Fij  to  Fmaxihj)- 

3.  Determine  the  minimum  gradient  magnitude  for  cells  near  highly  reacted  cells,  2  ■ 

a.  Calculate  the  lengths  l2,adj  and  h,diag- 

b.  Assign  the  characteristic  length  h  &om  h,adj  or  l^^iHag  depending  on  whether  Fmax.adjiiij) 
or  F„ax,diagiid)  is  greater  than  Ferit^ 

c.  Insure  that  the  result  of  the  calculation  is  not  less  than  the  smallest  spacing  between 
cell  (t,  j)  and  its  adjacent  neighbors. 

d.  Calculate  the  value  of  3 

6.  Given  |Vc|^,  i ,  |Vc|g,.^j  j  ,  and  modify  |Vc| ,  if  necessary. 

a.  If  either  F„ax  <  Fcru  and  jVcl  >  |Vc|^^,  i  or  F„ax  >  Ferit  and  |Vc|  >  |Vc|^^j  j  ’ 

|Vc|  remains  unchanged. 

b.  If  Fmax  <  Farit  and  |Vc|  <  2 ,  then  |Vc|  is  set  equal  to  zero. 

c.  If  Fnxx  >  F„n  and  lVc|  <  |Vc|,^<t  3  ,  then  |Vc|  is  set  equal  to  |Vc|^„t  3  . 

7.  Scale  |Vc|  back  to  the  local  temperature. 

Determination  of  model  constants 

The  gradient  model  contains  three  numerical  constants:  Ct,C2,  and  Ferit-  To  establish  the 
sensitivity  of  the  model  to  values  of  these  constants,  a  test  problem  was  formulated  on  a  uniform 
10  X  10  planar  grid.  The  cell  size  is  1.0  cm  square.  All  of  the  ceils  except  for  cell  (1,1)  contain  a 

stoichiometric  mixture  of  fuel  and  oxidant.  Cell  (1,1)  contains  the  amount  of  product  that  would 

result  from  complete  reaction  of  the  initial  mixture  in  any  of  the  other  cells.  This  simulates  a 
reaction  kernel  in  a  totally  premixed  reactant  field.  Since  there  is  no  convection  or  other  timestep 
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limiting  process  in  this  test,  the  numerical  timestep  may  be  chosen  arbitrarily  and.  as  such,  is  set 
to  1.0  s.  The  input  burning  velocity  is  chosen  to  be  0.01  cm/s,  so  that  the  time  for  the  front  to 
traverse  a  cell  is  long  enough  to  provide  adequate  time  resolution  of  the  test  results.  In  time,  the 
reaction  front  expands  and  the  radius  of  the  front  increases  by  one  cell  every  100  timesteps.  A  test 
problem  with  timestep  0.001  s  and  burning  velocity  10  cm/s  would  produce  identical  results.  Each 
of  the  tests  was  run  for  1000  steps. 

Before  proceeding  with  the  test  problem,  it  is  useful  to  define  some  terminology.  In  general, 
the  reaction  front  can  advance  across  a  cell  with  any  orientation  with  respect  to  the  planar  axes. 
Two  specific  orientations  are  at  the  focus  of  this  study:  horizontal  and  diagonal.  A  horizontal  cell 
is  one  for  which  the  reaction  front  advances  horizontally  across  the  cell.  For  the  test  problem  this 
implies  cells  j  =  1, . . . ,  10.  A  diagonal  cell  is  a  cell  in  which  the  front  crosses  diagonally,  i.e., 

cells  i  =  1,...,10.  Note  that  by  symmetry  vertical  cells  are  identical  to  the  horizontal  ones 

in  this  test  problem. 

Several  different,  specific  timesteps  in  a  calculation  axe  also  helpful  in  interpreting  the  results. 
That  timestep  during  which  a  cell  first  reacts  is  denoted  fr  and  that  timestep  when  reaction  is 
first  completed  in  a  cell  is  cr.  Tables  la  and  Ib  contains  /r  and  er  values  from  a  test  in  which 
Cl  =  C2  =  0.67  and  Fcru  =  0.95.  Note  that  a  given  cell  can  begin  reacting  before  reaction  is 
completed  in  its  neighbors.  Thus  the  reaction  region  will  extend  over  several  cells  at  a  given  time. 
The  timestep  fr  is  an  indicator  of  the  leading  edge  of  the  reaction  region  and  cr  in  an  indicator  of 
the  trailing  edge.  The  difference  between  fr  amd  cr  for  a  given  cell  is  the  cell  reaction  time.  (Since 
the  numerical  timestep  is  set  at  unity,  the  difference  between  fr  and  cr  is  not  only  the  number  of 
timesteps  the  calculation  has  proceeded  but  also  the  time  elapsed.)  The  steady-state  thickness  of 
the  reaction  region  is  not  well  established  until  the  horizontal  and  diagonal  cell  reaction  times  have 
reached  a  constant  value  along  their  respective  directions  of  the  front  advancement.  This  typically 
does  not  occur  until  reaction  is  complete  in  ail  the  cells  diagonally  adj2u:ent  to  the  reaction  kernel. 
After  the  reaction  region  has  been  established,  the  cell  reaction  times  are  approximately  equal  to 
the  ratio  of  the  thickness  of  the  reaction  region  and  the  input  flame  velocity.  For  the  example  in 
Tables  la  and  Ib,  the  horizontal  and  diagonal  cell  reaction  times  are  148.  Thus,  the  reaction  zone 
is  approximately  1.48  cells  thick.  This  test  is  further  discussed  below. 

Two  final  definitions  are  a  measure  of  the  burning  velocity  derived  from  the  calculations.  These 
quantities  provide  a  test  of  consistency  between  input  and  output.  The  difference  between  fr  (or 
cr)  for  one  horizontal  cell  and  fr  (or  cr)  for  its  adjacent  neighbor  cell  is  the  number  of  timesteps  it 
takes  the  reaction  front  to  advance  through  a  horizontal  cell.  This  quantity  is  called  the  horizontal 


advance  time.  The  counterpart  for  diagonal  cells,  the  diagonal  advance  time,  is  the  analagous 
difference  between  neighbors  which  are  diagonal  cells.  For  the  present  problem,  a  burning  velocity 
of  0.01  cm/s  and  a  cell  length  of  1.0  cm  implies  that  the  horizontal  advance  time  should  equal  100 
2uid  diagonal  advance  time  should  approximate  100  x  or  142. 

The  gradient  model  does  not  explicitly  define  the  location  of  the  reaction  front,  but  we  need 
such  a  definition  to  visualize  the  progress  of  a  reaction.  An  appropriate  tool  for  the  present  study 
is  a  contour  plot  of  the  er  values.  The  contours  generated  aire  called  reaction  completion  contours, 
and  the  location  of  the  front  at  any  given  timestep  is  defined  as  the  reaction  completion  contour 
generated  for  that  timestep.  It  should  be  emphasized  that  this  definition  of  a  front  location  is  not 
necessary  for  use  of  the  gradient  model  and  it  in  no  way  affects  the  results  of  the  model. 

To  begin  the  investigation  of  values  for  the  constants,  the  value  of  Feru  '^as  set  to  0.95  and  it 
was  assumed  that  Ci  =  =  C.  Typical  reaction  completion  contours  for  a  range  of  C  values  are 

shown  in  Figure  3.  Tables  of  fr  and  cr  such  as  those  shown  in  Tables  la  and  Ib  were  also  generated. 
A  summary  of  relevant  information  obtained  from  such  tables  for  different  values  of  C  is  shown  in 
Table  II. 

The  first  question  is  how  well  do  the  various  values  of  C  reproduce  the  true  horizontal  and 
diagonal  advance  times.  As  hoped,  the  horizontal  advance  time  is  independent  of  C  amd  the 
diagonal  advance  time  is  also  little  affected  except  for  C  >  0.7.  Thus  the  input  burning  velocity  is 
well  reproduced  over  a  large  range  of  C  values.  We  are  free  to  select  C  based  on  more  subtle  reasons 
than  simply  having  to  reproduce  the  desired  burning  velocity.  Table  II  shows  that  the  horizontal 
ceil  reaction  time  is  approximately  equal  to  1/(C  x  V&).  This  is  true  for  all  values  of  C.  The  diagonal 
cell  reaction  time  also  varies  with  1/(C  x  T4),  but  the  relationship  is  not  as  strong.  The  net  result  is 
that  the  curves  of  the  horizontal  and  diagonal  cell  reaction  times  cross  at  approximately  C  =  0.67. 

The  equality  of  the  horizontal  and  di^onal  cell  reaction  times  has  two  desirable  properties. 
First,  since  the  cell  reaction  time  is  proportional  to  the  reaction  region  thickness,  equality  of  the 
horizontal  and  diagonal  times  insures  that  the  reaction  region  has  a  uniform  thickness,  irregardless 
of  the  orientation  of  the  reaction  region.  Second,  once  the  reaction  region  has  become  established, 
equality  of  the  times  insures  the  the  reaction  region  will  propagate  with  little  additional  distortion. 

The  contour  plots  in  Figure  3  reinforce  the  choice  of  C  =  0.67.  Recall  that  for  the  test  problem 
the  reaction  front  should  propagate  as  concentric  circles  after  the  initial  square  cell  has  smoothed 
out.  Even  though  the  horizontal  and  diagonal  advance  times  equal  their  true  values  for  C  <  O.S,  the 
contours  a^e  most  accurately  circular  for  C  near  0.67.  The  cause  of  the  flattening  for  C  <  0.67  can 
be  traced  to  the  time  needed  to  completely  react  the  horizontal  and  diagonal  cells  directly  connected 


to  the  reaction  kernel  cell.  The  values  of  cr  for  cells  (1,2)  and  (2,2)  are  also  shown  in  Table  II  along 
with  the  difference  between  them.  For  the  smaller  values  of  C,  the  difference  indicates  that  the 
front  initially  propagates  much  more  quickly  in  the  horizontal  and  vertical  directions  than  along 
the  diagonals.  The  result  is  the  perceived  flattening.  Because  the  diagonal  cell  reaction  times  are 
less  than  their  horizontal  counterparts  for  these  values  of  C,  the  reaction  fronts  would  eventually 
appear  more  circular  after  an  additional,  significant  delay. 

The  initial  test  cases  assumed  Ci  =  Cj.  The  next  set  of  tests,  shown  in  Tables  III  and  IV, 
assume  Ci  #  Cj.  Again  Fcru  is  set  to  0.95.  The  most  obvious  result  is  that  the  horizontal  advance 
time  equals  its  true  value  only  if  C\  and  are  equal.  Thus,  the  initial  assumption  that  C\  =  Ci 
appears  to  be  an  optimal  condition  for  the  model. 

Finally,  the  effects  of  varying  Ferit  were  investigated,  as  shown  in  Table  V.  Here  Ci  =  C2  = 
0.67.  The  result  is  that  Ferit  =  0.95  appears  to  be  optimum  although  distortions  introduced  by  a 
slightly  different  value  are  minor.  Varying  Ci  and  Ci  to  improve  the  results  for  a  different  value  of 
Ferit  did  not  produce  an  optimum  C  value  appreciably  different  from  0.67. 

Selection  of  burning  velocity 

The  burning  velocity  is  an  input  to  this  model.  It  may  be  specified  as  a  constant  value  or  may 
be  a  function  of  any  relevant  parameters.  For  e.xample,  could  be  function  of  cell  composition  or 
temperature.  The  choice  of  burning  velocity  is  totally  independent  of  the  front  propagation  model. 

Change  of  species  and  energy  release  due  to  reaction 

Given  |Vc| ,  the  change  in  species  number  density  amd  the  energy  release  can  be  calculated. 
The  change  in  the  number  density  of  the  product  for  a  cell  is 

Ac  =  VfcAf|Vc|.  (24) 

It  is  important  to  insure  that  more  product  is  not  created  than  there  are  reactants  available.  The 
change  in  fuel  number  density  is  first  calculated  as  the  minimum  of  the  amount  of  fuel  needed  to 
produce  the  change  in  number  density  of  the  product  or  the  amount  of  fuel  present.  The  same  is 
true  for  the  oxidant.  Then 


where  /  and  o  aire  the  fuel  and  oxidant  number  densities,  respectively,  iV/  and  No  are  the  number 
of  fuel  and  oxidant  molecules  consumed  per  reaction,  and  Nc  is  the  number  of  product  molecules 
produced  per  reaction.  The  actual  species  changes  in  number  density  are  computed  to  assure 
consistency  and  mass  conservation, 


(A/y  =  min 
(Ao)'  =  min 
(Ac/  =  min 
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(26a) 

(266) 

(26c) 


The  energy  change  is 

This  procedure  guarantees  mass  conservation  and  appropriate  energy  release  due  to  reaction. 


A  (^«)' 

Ae  =:  ^-rj: —  x  energy  per  reaction. 

iVg 


(27) 


Testing  the  Gradient  Model 


An  extensive  series  of  tests  were  performed  to  see  if  the  gradient  method  reproduces  the 
properties  of  a  reaction  front  that  propagates  at  a  known  velocity.  This  is  a  consistency  test 
between  input  and  output.  For  reacting  but  nonconvecting  media,  the  reaction  completion  contours 
introduced  above  au:e  used  to  track  the  progress  of  the  front.  For  reactions  occurring  concurrently 
with  convection,  the  reaction  front  is  defined  as  the  0.80  contour  of  the  Fij  values.  All  of  the 
following  tests  were  performed  with  Ci  =  C2  =  0.67  and  Feru  =  0.95. 

The  first  test  was  an  extension  of  the  test  used  to  determine  the  model  constants.  Here  a 
20  X  20  uniform  grid  was  initialised  with  a  stoichiometric  mixture  of  fuel  amd  oxidant  in  all  cells 
except  for  cell  (10,10).  Cell  (10,10)  contains  the  amount  of  product  that  would  result  from  complete 
reaction  in  any  of  the  other  cells.  There  is  no  convection  present  in  this  test.  The  test  was  run 
for  1500  steps  and  the  results  are  shown  in  Figures  4  and  5.  The  reaction  completion  contours  in 
Figure  4  show  that  the  reaction  front  propagates  away  from  the  reaction  kernel  in  concentric  circles 
and  the  spacing  between  the  circles  indicates  that  the  front  is  propagating  at  the  correct  speed. 
The  initial  raggedness  in  the  contours  is  due  to  the  coarseness  of  the  test  grid.  The  innermost 
contour  is  diamond-shaped  because  the  initial  condition  for  the  test  consists  of  a  step  discontinuity 
at  the  interface,  and  the  gradient  method  must  propagate  the  interface  through  the  cells  adjacent 
to  the  initial  reaction  kernel  before  the  interface  profile  can  be  established.  The  first  contour  for 
which  the  interface  profile  is  established  is  the  contour  for  step  200.  It  should  be  emphasized  that 
the  finite  thickness  of  the  profile  is  the  numerical  limit  of  resolution  and  is  not  a  profile  connected 
with  the  physical  interface. 

Figure  5  shows  the  time  the  front  reaches  a  given  distance  from  the  reaction  kernel.  Time  on 
the  vertical  axis  is  measured  by  cr  and  distance  is  measured  from  the  center  of  the  reaction  kernel 
to  the  cell  vertex  farthest  from  the  reaction  kernel.  This  distance  is  chosen  because  reaction  is  not 
complete  until  the  physical  interface  has  reached  this  farthest  point.  The  solid  diagonal  line  in 
Figure  5  is  the  true  position  of  an  interface  initiated  at  a  point  source  located  at  the  center  of  the 
actual  reaction  kernel  and  propagating  with  the  speed  V^.  The  data  points  are  the  results  of  the 
test  calculation.  There  are  two  reasons  for  the  small  displacement  between  the  solid  line  and  the 
data.  First,  the  actual  reaction  kernel  is  not  a  point  source  and  interface  propagation  begins  on 
the  square  boundary  of  the  reaction  kernel.  This  explains  the  data  point  on  the  x-axis  displaced 
from  the  origin,  and  that  displacement  is  an  upper  estimate  of  the  horizontal  displacement  of  the 
other  points.  Adjustment  for  the  finite  reaction  kernel  moves  the  data  points  to  the  right.  Second, 
there  is  a  delay  in  the  reaction  completion  of  the  cells  adjacent  to  the  reaction  kernel.  This  delay  is 


related  to  the  time  it  takes  to  establish  the  steady-state  interface  profile,  and  so  all  completion  times 
are  slightly  longer  than  if  the  interface  started  with  its  steady-state  propagation  rate.  Correction 
for  this  effect  moves  the  data  points  down.  The  exact  magnitude  of  each  correction  depends  on 
the  location  of  the  point  of  interest  with  respect  to  the  reaction  kernel.  If  all  points  are  corrected 
by  the  displacement  of  the  data  on  the  x-axis  and  the  time  delay  for  reaction  completion  of  the 
cells  hori2ontally  adjacent  to  the  reaction  kernel  cell,  then  the  data  points  fall  almost  exactly  on 
the  point  source  solution  line.  These  finite  resolution  corrections  do  not  alter  the  slope  of  a  line 
suggested  by  the  data  points  in  Figure  5,  and  this  slope  is  almost  identical  to  the  input  propagation 
speed. 

The  second  test  investigated  the  ability  of  the  gradient  model  to  handle  merging  fronts.  A 
40  X  40  uniform  grid  was  used  and  two  reaction  kernels  were  specified.  The  test  was  run  for  2100 
steps.  Figure  6  shows  that  the  fronts  merge  without  difficulty.  The  reaction  fronts  propagate 
independently  until  they  merge,  and  then  the  merged  front  continues  to  propagate  as  the  union 
of  two  independent  fronts.  The  cusps  seen  at  the  points  where  the  fronts  merge  occur  because  a 
constant  propagation  speed  was  specifed  and  no  attempt  was  made  here  to  model  the  finer  details 
of  the  reaction.  Several  combinations  of  point  and  line  reaction  kernels  have  been  tested  in  both 
flowing  and  stationary  backgrounds,  and  the  gradient  model  always  handles  interface  merging  well. 

The  next  test  case  exercised  the  gradient  model  on  a  nonuniform  grid.  The  grid  is  uniformly 
stretched  left-to-right  (x-direction)  by  a  factor  of  1.03  and  uniformly  stretched  bottom-to-top  (y- 
direction)  by  a  factor  of  1.07.  The  crosshatched  rectangles  in  the  comers  of  Figure  7  show  the  cell 
sizes  and  shapes  at  those  locations.  For  this  test,  a  single-cell  reaction  kernel  was  initialized  in  a 
premixed  reactant  field.  The  test  was  run  for  2400  steps.  Reaction  completion  contours  hte  shown 
in  Figure  7.  The  distortion  in  the  three  innermost  contours  occurs  because  the  interface  profile 
is  not  yet  established.  Any  remaining  deviations  from  a  circular  trace  relax  as  the  calculation 
proceeds.  A  graph  of  time  as  a  function  of  distance  from  the  reaction  kernel  is  very  similar  to 
Figure  5,  and  again  the  slope  suggested  by  the  data  is  nearly  identical  to  the  specified  propagation 
speed. 

The  remaining  tests  were  performed  in  convecting  flows.  Since  the  gradient  model  keys  on 
sharp  gradients  in  the  flow,  it  is  important  that  the  convection  algorithm  is  not  highly  diffusive  at 
the  front,  and  any  ‘*wiggles’’  that  distort  the  local  gradient  should  not  appear  in  the  solution  of  the 
indicator  variable.  Here  the  Flux- Corrected  Transport  (FCT)  algorithm  [11]  is  used  to  transport 
mass,  momentum,  energy,  and  species  in  fast  flows,  and  BIC-FCT  [12],  an  implicit  corrector  used 
with  the  FCT  algorithm,  is  used  to  convect  slower  flows.  The  expansion  due  to  energy  releeise 


follows  naturally  from  a  source  term  in  the  energy  conservation  equation  and  does  not  require  the 
calculation  of  a  separate  velocity  field  as  is  done,  for  example,  in  Ref.  13.  Examples  using  both 
FCT  algorithms  will  now  be  presented. 

The  first  test  with  convection  was  performed  on  a  uniform  grid  with  an  x-component  velocity 
of  5  m/s  and  a  y-component  velocity  of  10  m/s.  A  3  x  3  reaction  kernel  was  initialized  in  a  premixed 
reactant  field.  The  burning  velocity  was  specified  as  10  m/s.  The  standcird  FCT  algorithm  was 
used  to  simulate  convection,  A  series  of  Fij  =  0.80  contours  are  shown  in  Figure  8.  The  innermost 
contour  represents  the  reaction  front  at  timestep  100  and  each  successive  contour  shows  the  reaction 
front  100  timesteps  later.  All  sides  of  the  reaction  front  move  with  speeds  which  are  the  vector  sum 
of  the  background  velocity  and  the  propagation  speed  in  the  direction  of  the  outward  normal.  This 
results  in  the  lower  edge  of  the  interface  remaining  stationary  because  the  downward  propagation 
just  balances  the  fiow.  Note  also  that  the  initially  square  reaction  front  is  becoming  progressively 
more  circular.  Tests  with  convection  and  several  merging  interfaces  show  the  same  front  behavior. 

All  of  the  previous  tests  were  done  without  energy  release.  The  final  convection  test  involves 
a  comparison  of  the  evolution  of  the  flow  caused  by  a  single  vortex  in  a  box  when  there  is  no 
reaction,  when  the  reaction  is  isothermal,  and  when  there  is  reaction  with  energy  release.  The  flow 
is  initialized  with  the  velocities  in  the  central  “box"  of  a  3  x  3  array  of  potential  point  vortices. 
By  locating  the  center  of  the  central  vortex  at  the  intersection  of  computational  grid  lines,  the 
calculation  of  velocities  at  the  cell  centers  avoids  the  singularity  at  the  vortex  core.  The  circulation 
of  each  vortex  is  10000  cm^/s.  A  40  cm  square  computational  domain  is  used  and  the  grid  is 
slightly  stretched  to  concentrate  cells  In  the  center  of  the  box.  The  burning  velocity  is  200  cm/s 
and  the  timestep  is  3.75  x  10“®  s.  The  upper  half  plane  of  the  box  is  initialized  with  a  stoichiometric 
mixture  of  fuel  and  oxidant  and  the  lower  half  plane  contains  only  product.  The  temperature  of 
the  entire  field  is  1000®K,  except  for  the  energy  release  case  where  the  temperature  of  the  products 
is  1950^K.  The  pressure  is  uniform  at  I  atm.  BIC-FCT  was  used  for  convection. 

Figure  9  shows  contours  of  Fij  at  intervals  of  0.1.  For  the  nonreacting  case,  the  vortex 
grows  with  time,  but  the  detail  of  the  vortex  center  is  smeared  due  to  numerical  diffusion.  This 
smearing  is  the  result  of  the  coarseness  of  the  test  grid.  Otherwise,  the  velocity  near  the  vortex 
center  is  sufficiently  high  to  produce  several  windings  of  the  central  spiral.  However,  even  with 
the  numerical  smearing,  the  interface  is  still  symmetrically  located  about  the  vortex  center.  In 
the  isothermal  reaction  case,  several  differences  are  seen.  First,  the  propagation  model  tends  to 
sharpen  the  interface.  The  distance  between  the  0.1  and  0.9  contours  at  the  edges  of  the  box 
has  decreased  because  the  propagation  algorithm  has  reacted  numerically  diffused  material.  This 


effect  is  most  obvious  near  the  vortex  center  where  the  reaction  front  is  progressing  into  the  vortex 
because  the  rate  of  reaction  is  exceeding  the  rate  of  numerical  diffusion.  The  contours  are  not 
symmetric  because  1)  the  reaction  model  is  only  active  on  the  product  side  of  the  interface  and  so 
interface  sharpening  cannot  occur  on  the  reactants  side,  and  2)  the  reaction  front  is  propagating 
into  the  reactants  arm  of  the  vortex  spiral.  Because  there  is  no  energy  release  in  the  isothermal  case, 
the  mass,  momentum,  and  energy  is  unchanged  from  the  nonreacting  case.  Only  the  individual 
species  concentrations  are  altered  by  the  reaction.  This  is  not  the  situation  when  the  material 
conversion  is  accompanied  by  energy  release.  In  that  case  the  vortex  flow  is  augmented  due  to  the 
expansion,  and  this  causes  the  reactants  not  only  to  be  consumed  due  to  reaction  but  also  pushed 
away  from  the  neighborhood  of  the  vortex  center  into  lower  velocity  parts  of  the  box.  The  interface 
sharpening  effect  is  still  present.  Together,  the  continued  presence  of  the  interface  sharpening  and 
the  expansion  in  a  direction  away  from  the  reaction  location  indicate  that  the  energy  deposition 
is  being  done  correctly.  In  all  convecting  cases,  the  interface  moves  with  the  sum  of  the  flow  and 
propagation  velocities. 
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Discussion 


The  gradient  method  has  been  developed  to  propagate  an  interface  which  is  characterized 
by  large  gradients  due  to  conversion  of  material  at  the  interface.  Although  the  development  and 
testing  described  above  were  done  for  two-dimensional  problems,  the  formulation  can  easily  be 
extended  to  three  dimensions.  Unlike  other  interface  tracking  methods,  the  gradient  method  does 
not  require  the  location  and  movement  of  the  interface  to  be  defined  e.xplicitly.  Only  the  effects  of 
the  interface  are  considered.  The  gradient  method  has  demonstrated  its  ability  to  simulate  reaction 
at  an  interface  on  uniform  or  stretched  grids.  There  is  no  limit  to  the  number  of  interfaces  which 
can  be  present  and  any  number  of  merging  interfaces  can  be  handled  without  difficulty. 

This  version  of  the  gradient  method  uses  the  product  number  density  as  the  indicator  variable, 
bat  other  versions  are  possible  using  other  indicators.  For  example,  assuming  that  expansion  due  to 
energy  release  causes  proportionate  amounts  of  the  species  in  a  cell  to  leave  that  cell  indicates  that 
using  the  product  volume  fraction  as  the  indicator  would  eliminate  the  need  of  scaling  to  a  reference 
temperature.  However,  this  choice  may  not  reduce  the  total  number  of  arithmetic  operations  used. 
The  approach  presented  in  this  paper  is  sufficiently  general  to  accommodate  a  variety  of  alternate 
indicator  variables. 

As  well  as  the  gradient  method  has  been  shown  to  perform,  there  are  still  some  limitations, 
many  of  which  also  apply  to  the  volume  tracking  methods.  Because  the  location  of  the  interface 
is  known  only  approximately,  information  such  as  the  curvature  of  the  interface  must  be  derived 
from  a  separate  calculation.  For  example,  Ghoniem  [14]  has  recently  used  an  approach  similar 
to  Chorin’s  idea  of  osculating  circles  [15]  to  estimate  the  curvature  of  a  reaction  front  defined  by 
SLIC  [4,16].  A  more  important  limitation  of  any  method  that  reconstructs  the  interface  on  the 
basis  of  a  marker  or  indicator  quantity  is  that  care  must  be  taken  to  insure  that  the  presence  of 
the  indicator  in  a  cell  is  solely  the  result  of  reaction  in  that  cell.  For  example,  small  amounts  of 
the  indicator  numerically  diffused  from  the  reaction  front  should  not  result  in  the  identification  of 
spurious  interfaces.  This  care  is  necessary  regardless  of  the  choice  of  marker  or  indicator  variable. 

The  concern  over  the  extraneous  appearance  of  the  indicator  away  from  the  interface  highlights 
the  need  to  choose  a  convection  algorithm  which  avoids  excessive  numerical  diffusion.  Additionally, 
the  convection  algorithm  should  be  monotonic  so  that  no  wiggles  appear  in  the  indicator  profile 
which  will  introduce  large  errors  in  the  estimation  of  the  gradient.  The  present  implementation  has 
tested  two  variations  of  the  FCT  algorithm  to  satisfy  this  requirement. 

Finally,  an  important  requirement  for  any  algorithm  that  is  to  be  used  in  a  numerical  simulation 
is  that  it  allows  full  use  of  vector  and  parallel  computing  systems.  The  program  logic  of  the  new 
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method  vectorizes  fully.  In  addition,  if  the  chosen  indicator  variable  is  already  included  in  the 
calculation,  then  no  extra  array  storage  is  necessary  beyond  whatever  scratch  space  is  useful  to 

enhance  vectorization. 
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Table  la.  TImestep  of  first  reaction,  fr. 


J\I 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 

0 

1 

92 

190 

289 

389 

488 

588 

688 

788 

2 

1 

1 

119 

212 

308 

405 

502 

600 

699 

798 

3 

92 

119 

150 

250 

345 

437 

531 

626 

722 

819 

4 

190 

212 

250 

305 

387 

481 

571 

663 

756 

850 

5 

289 

308 

345 

387 

452 

528 

618 

707 

797 

888 

6 

389 

405 

437 

481 

528 

596 

670 

755 

844 

933 

7 

488 

502 

531 

571 

618 

670 

738 

812 

894 

982 

8 

588 

600 

626 

663 

707 

755 

812 

881 

954 

*** 

9 

688 

699 

722 

756 

797 

844 

894 

954 

*** 

10 

788 

798 

819 

850 

888 

933 

982 

*** 

*** 

Table  Ib. 

Timestep  when  reaction 

is  completed, 

cr. 

J\I 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 

0 

137 

237 

337 

436 

536 

636 

736 

836 

936 

2 

137 

154 

255 

353 

451 

549 

647 

746 

845 

946 

3 

237 

255 

309 

393 

486 

580 

675 

771 

867 

966 

4 

337 

353 

393 

458 

534 

623 

714 

807 

900 

997 

5 

436 

451 

486 

534 

601 

675 

761 

851 

941 

6 

536 

549 

580 

623 

675 

744 

818 

900 

988 

7 

636 

647 

675 

714 

761 

818 

886 

960 

*** 

*** 

8 

736 

746 

771 

807 

851 

900 

969 

*** 

*** 

*** 

9 

836 

845 

867 

900 

941 

988 

*** 

*** 

10 

936 

946 

966 

997 

««* 

*** 

*** 

*** 

*** 

0** 

Table  n.  Derived  values  for  Ci  =  C2  =  C  and  Ferit  —  0.95. 


c 

horizontal 

diagonal 

horizontal  cell 

diagonal  cell 

cr 

cr 

cr(2,2)- 

advance  time 

advance  time 

reaction  time 

reaction  time 

(1,2) 

(2,2) 

cr(l,2) 

0.40 

100 

143 

248 

222 

153 

215 

62 

0.50 

100 

143 

199 

185 

149 

190 

41 

0.60 

100 

142 

166 

153 

144 

166 

22 

0.87 

100 

142 

148 

148 

137 

154 

17 

0.70 

100 

142 

141 

147 

133 

150 

17 

0.80 

100 

132 

124 

138 

122 

138 

16 

0.90 

100 

122 

111 

127 

111 

131 

20 

1.00 

96 

111 

100 

115 

101 

123 

22 

21 


PBPU  WJ' 


Table  III.  Derived  values  for  Ci  variable,  C2  =  0.67,  and  Fcrxt  =  0.95. 


Cl 

horizontal 

diagonal 

horizontal  cell 

diagonal  cell 

CT 

CT 

cr(2,2)- 

advance  time 

advance  time 

reaction  time 

reaction  time 

(1,2) 

(2,2) 

cr(  1.2) 

0.50 

88 

133 

161 

163 

140 

165 

25 

0.60 

94 

142 

152 

153 

140 

158 

18 

0.67 

100 

142 

148 

148 

137 

154 

17 

0.70 

104 

143 

144 

148 

135 

153 

18 

0.80 

111 

143 

135 

148 

131 

152 

21 

Table  IV.  Derived  values  for  Ci  =  0.67,  Cj  vnriAble,  and  Fcrit  =  0.95. 


c. 

horizontal 

diagonal 

horizontal  cell 

diagonal  cell 

cr 

CT 

cr(2.2)- 

advance  time 

advance  time 

reaction  time 

reaction  time  (1,2) 

(2,2) 

cr(l,2) 

0.50 

107 

142 

156 

148 

136 

178 

42 

0.60 

105 

143 

153 

148 

139 

162 

23 

0.67 

100 

142 

148 

148 

137 

154 

17 

0.70 

97 

141 

143 

147 

135 

151 

16 

0.80 

89 

130 

129 

136 

123 

141 
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Table  V.  Derived  values  for  Ci  =  C2  =  0.67  and  fcrit  =  F. 


F 

horizontal 

diagonal 

horizontal  cell 

diagonal  cell 

cr 

cr 

cr(2,2)- 

advance  time 

advance  time 

reaction  time 

reaction  time 

(1,2) 

(2,2) 

cr(l,2) 

0.80 

100 

132 

148 

160 

137 

157 

20 

0.90 

100 

141 

148 

154 

137 

154 

17 

0.95 

100 

142 

148 

148 

137 

154 

17 

0.99 

100 

142 

148 

143 

137 

154 

17 

22 


Figure  6.  Reaction  completioa  contours  for  merging  interfaces. 
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Figure  8.  Fi,i  =  0.80  contours  (or  convection  test. 


Figure  9.  Fi^j  contours  for  vortex  tests.  (Topmost  contour,  Ffj  =  0.1;  bottommost  contour, 


